Interactive Deforestation Map

This document can be used to visualize a single protected area from the World Database on Protected Areas (WDPA) on an interactive web map. The map displays the selected protected areas (PAs) including a buffer region around it. The barplot shows the amount of forest cover loss over time for the PA itself as well as for the buffer zone around it. You can adjust the params object in the YAML header to your requirements.

Code
basemap_custom <-
  leaflet() |>
  # add external map providers
  addTiles(group = "OpenStreetMap") |>
  addProviderTiles(providers$Esri.WorldImagery, group="Satellite") |>
  addProviderTiles(providers$CartoDB.Positron, group="CartoDB") |>
  addProviderTiles(providers$Esri.WorldShadedRelief, group="Topography") |>
  addProviderTiles(providers$NASAGIBS.ViirsEarthAtNight2012, group="Nightlights") |>
  addTiles(
    "https://tiles.globalforestwatch.org/umd_tree_cover_loss/latest/dynamic/{z}/{x}/{y}.png",
    group = "Forest Cover Loss (2001-2020)",
    attribution = "Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. “High-Resolution Global Maps of 21st-Century Forest Cover Change.” Science 342 (15 November): 850–53. Data available on-line from: http://earthenginepartners.appspot.com/science-2013-global-forest."
  ) |>
  addTiles(
    "https://tiles.globalforestwatch.org/umd_regional_primary_forest_2001/latest/dynamic/{z}/{x}/{y}.png",
    group = "Regional Primary Forests (2001)",
    attribution = "Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. “High-Resolution Global Maps of 21st-Century Forest Cover Change.” Science 342 (15 November): 850–53. Data available on-line from: http://earthenginepartners.appspot.com/science-2013-global-forest."
  ) |>
  addFullscreenControl() |>
  # add legend(s)
  addLayersControl(
    baseGroups = c("Satellite", "CartoDB", "OpenStreetMap", "Topography", "Nightlights"),
    overlayGroups = c("Protected Area", "Buffer Zone", "Forest Cover Loss (2001-2023)", "Regional Primary Forests (2001)"),
    options = layersControlOptions(collapsed = FALSE)) |>
  # uncheck some layers in layer control
  hideGroup(group = c("Regional Primary Forests (2001)","Labels (PA Names)"))

deforestation_map <- basemap_custom |>
  addPolygons(data = buffer, fillOpacity = 0.0, opacity = 0.8, color = params$col_bf, smoothFactor = 0, weight =params$lwd_bf, group = "Buffer Zone",
              dashArray = "10 8") |>
  addPolygons(data = wdpa, fillOpacity = 0.0, opacity = 1.0, color = params$col_pa, smoothFactor = 0, weight = params$lwd_pa, group = "Protected Area")

Loss of forest cover over time

Code
plt <- ggplot(data = inds) +
  geom_col(aes(x=datetime, y=losses, fill=buffer), position=ifelse(params$stacked, "stack", "dodge")) +
  labs(x = "Year", y = "Forest cover losses (in ha)", fill = "Zone") +
  scale_x_datetime(date_breaks = "1 year", date_labels =  "%Y") +
  scale_fill_manual(values = c(params$col_pa, params$col_bf)) +
  theme_classic() +
  theme(
    axis.title=element_text(size=16),
    axis.text=element_text(size=14),
    axis.text.x=element_text(angle=60, hjust=1),
    legend.position="bottom"
    ) 

ggplotly(plt)